function [M,PI]=mdrkw_learning(B,CF,Q,RHO,LE,MU,reord,nd,lpd,nx,bcrit,schur,bigA)

% This program takes the reduced system produced by redkw.m and 
% the decomposition of its dynamic subsystem undertaken by dynkw.m
% and computes the RE solution. Computation can take place using either
% the eigenvalue-eigenvector decomposition or the Schur decomposition

% relative to mdrkw, solves the problem when exogenous shocks are not
% perfectly observed but instead Kalmn filtering is used to learn about
% them (Kevin Moran, 2012-2015)


v=version;

if (eval(v(1))>=5)
   npd=length(lpd);
else
   npd=max(size(lpd));   % npd is number of predetermined variables
end

nus=rows(LE); % number of unstable roots.


[ind back]=sort(reord); % reverse transformation

%if (nargin<13)
   toli=10^-6;
%end

if (nus>0)
   
ny=rows(B);
nf=ny-nd;


ndrv=cols(Q);

% Consider the equation: E d(t+1|t) = W d(t) + PSId(F) E x(t|t)

% The solution strategy is then as follows:  We know that 
%  (1)  LE*Ed(t+1|t) = LE*W d(t) + LE*PSId(F) Ex(t|t)
% We know further that LE*W = MU * LE from the definition
% of left eigenvectors (this holds for any set of eigenvectors).
% Hence, equation (1) is the unstable canonical variables 
% equation.  We can hence solve it "forward" as
%  (2)  LE*d(t) = inv(MU)*LE Ed(t+1|t) -inv(MU)*LE*PSId(F) Ex(t|t)

nlead=cols(CF)/nx;
nlead=nlead-1;

% Calculate CF0*Q*I+CF1*Q*RHO+CF2*Q*RHO^2 + ...CFn*Q*RHO^nlead
% i.e., calculate C(F) Ex(t|t)
R=eye(cols(RHO));
PHI=CF(:,1:nx)*Q*R;
j=0;

add = PHI;
%see my notes for explanation (Kevin Moran)


while (j<nlead)
 j=j+1;
 R=R*RHO;
 PHI=PHI+CF(:,j*nx+1:(j+1)*nx)*Q*R;
end


% We now have a system in the form 
%   A* E y(t+1)|I(t) = B y(t) + PHI drv(t),
% with the A and B matrices having been "reduced" by dynkw.m
% but we still need to rule out explosive paths, etc. Part of this system
% is E d(t+1)|I(t) = W  d(t) +   PHId drv(t).  Multiplying this by LE we get
%    Eu(t+1)|I(t) =  MU u(t) + LE*PHId drv(t), which are "separated"
% dynamic equations for the nonpredetermined/unstable canonical
% variables u(t) since MU is upper triangular. Set q = LE*PHI2.

% To impose stability, we seek a matrix theta such that
% theta*RHO = MU * theta + q.  There are then two cases
% If MU is diagonal, then the ith row of theta satisfies:
% thetai * RHO = MU(i,i) * thetai + qi => thetai = qi/(RHO-MU(i,i)*eye(RHO);
% If MU is lower triangular, then we must add in some 
% additional terms that involve MU(i,j)*thetaj for j
% =1,2,...i-1.

q = LE*PHI(nf+1:ny,:);
qadd = LE*add(nf+1:ny,:);  %again see notes

bigB = eye(ndrv)-bigA;
% Look at notes for deficitions of "bigA" and "bigB" 
% (they arise from updating the Kalman filter)

theta1=zeros(nus,cols(PHI));
theta2=zeros(nus,cols(PHI)); % these go under the "nu" notation in the notes
IRHO=eye(cols(RHO));

for i=1:nus;
 theta1(i,:)=q(i,:);
 theta2(i,:)=q(i,:);
 if schur==1
   for j=1:i-1;
     theta1(i,:)=theta1(i,:)+MU(i,j)*theta1(j,:);
     theta2(i,:)=theta2(i,:)+MU(i,j)*theta2(j,:);
   end
 end
 theta1(i,:)= ( theta1(i,:)/(RHO-MU(i,i)*IRHO) )*bigA - (1/MU(i,i))*qadd(i,:)*bigB;
 theta2(i,:)= (theta2(i,:)/(RHO-MU(i,i)*IRHO)) * bigB + (1/MU(i,i))*qadd(i,:)*bigB;
end

% We now know that u(t) = LE * d(t) = theta * drv(t).  We solve for lam(t),
% the first nus elements of d(t), as lam(t) = -(Ll\Lk) k(t) + (Ll\theta) drv(t):

if (npd>0)
 lamk=-LE(:,1:nus)\LE(:,nus+1:nd);  
end

lamdr1= LE(:,1:nus)\theta1;
lamdr2= LE(:,1:nus)\theta2;



BT=B;

% We want to impose this solution on the dynamic system, i.e.,
% "substitute out" for lam(t) using the above rule, which we write
% as lam(t) - lamk * k(t) - lamdr * drv(t).  We have equations of the form:
%  B2 * lam(t) +         B3 * k(t) +   PHI           * drv(t) = 0.  Adding  
% -B2 * lam(t) +    B2*lamk * k(t) +   B2*lamdr      * drv(t) = 0,  we get
%  0  * lam(t) + B3+B2*lamk * k(t) +  (PHI+B2*lamdr) * drv(t) = 0

B2 = BT(:,nf+1:nf+nus);
B3 = BT(:,nf+nus+1:ny);
% incorporate the influence of k(t) via lam(t) into system.
if (npd>0)
 BT(:,nf+nus+1:ny)=B3+B2*lamk;
end
% incorporate the influence of drv(t) via lam(t) into system.
PHI1=PHI*bigA+B2*lamdr1+add*bigB;
PHI2=PHI*bigB+B2*lamdr2-add*bigB;

% patch up the lam(t) influences and the lam(t) equations.
BT(:,nf+1:nf+nus)=zeros(rows(B2),cols(B2));
BT(nf+1:nf+nus,nf+1:nf+nus)=eye(nus);
if (npd>0)
 BT(nf+1:nf+nus,nf+nus+1:nf+nd)=-lamk;
end

PHI1(nf+1:nf+nus,:) = -lamdr1;
PHI2(nf+1:nf+nus,:) = -lamdr2;

% solve for state space system:
% z(t) =   PI*s(t)
% s(t+1) = M s(t) + e(t+1)
% with s(t) = |k(t)  |   and z(t) = |y(t)|
%             |drv(t)|              |k(t)|
%                                   |x(t)|

ndrv=cols(Q);

PI= [-BT(1:nf+nus,nf+nus+1:ny) -PHI1(1:nf+nus,:) -PHI2(1:nf+nus,:)
      eye(npd)                  zeros(npd,ndrv)   zeros(npd,ndrv)  ];

% This is a PI matrix for the new ordering, we now create one
% for the old ordering and add exogenous variables. The index
% "back" reverses the earlier ordering.

PI=  [PI(back,:)
      zeros(nx,npd)              Q          zeros(nx,ndrv)  ];

M = [BT(nf+nus+1:ny,nf+nus+1:ny)  PHI1(nf+nus+1:ny,:)  PHI2(nf+nus+1:ny,:)
     zeros(ndrv,npd)              RHO                  zeros(ndrv,ndrv)
     zeros(ndrv,npd)              RHO*bigA             RHO*bigB            ];
 
 
%PI_l = PI;
%M_l = M;


   
   


else % there is a system with no unstable roots

ny=rows(B);
nf=ny-nd;

% Consider the equation: E d(t+1|t) = W d(t) + PSId(F) Ex(t|t)

% The solution strategy is then very simple because there
% are no variables to be solved "forward". W is simply the 
% feedback matrix and the PSI(F) polynomial governs
% the response to driving variables.


nlead=cols(CF)/nx;
nlead=nlead-1;

% Calculate CF0*Q*I+CF1*Q*RHO+CF2*Q*RHO^2 + ...CFn*Q*RHO^nlead
% i.e., calculate C(F) Ex(t|t)

R=eye(size(RHO));
PHI=CF(:,1:nx)*Q*R;
j=0;

while (j<nlead)
 j=j+1;
 R=R*RHO;
 PHI=PHI+CF(:,j*nx+1:(j+1)*nx)*Q*R;
end

% We now have a system in the form 
%   E y(t+1|t) = B y(t) + PHI drv(t)
% The first ny-np of these equations are for flows
% and are in the form I f(t) + V k(t) + PHIf drv(t)
% and the last np are in the form: k(t+1)=W k(t) +PHId drv(t)


% solve for state space system:
% x(t) =   PI*s(t)
% s(t+1) = M s(t) + e(t+1)
% with s(t) = |k(t)  |   and z(t) = |y(t)|
%             |drv(t)|              |k(t)|
%                                   |x(t)|

ndrv=cols(Q);

PI= [-B(1:nf,nf+1:ny) -PHI(1:nf,:)
      eye(npd)         zeros(npd,ndrv)];

% This is a PI matrix for the new ordering, we now create one
% for the old ordering and add exogenous variables. The index
% "back" reverses the earlier ordering.

PI=  [PI(back,:)
      zeros(nx,npd)              Q           ];

M = [B(nf+1:ny,nf+1:ny)   PHI(nf+1:ny,:)
     zeros(ndrv,npd)                 RHO];


end

% test for the presence of larger imaginary components

IPI=imag(PI);
t=max(max(abs(IPI)));

if t>toli;
   disp(['imaginary component of PI is larger than ', num2str(toli)])
   disp('This tolerance can be adjusted by setting the last argument')
   disp('of the mdrkw.m function, as called in resolkw.m')
   disp('   strike any key to continue')
   pause
end

IM=imag(M);
t=max(max(abs(IM)));

if t>toli;
   disp(['imaginary component of M is larger than ', num2str(toli)])
   disp('This tolerance can be adjusted by setting the last argument')
   disp('of the mdrkw.m function, as called in resolkw.m')
   disp('   strike any key to continue')
   pause
end
M=real(M);
PI=real(PI);

